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chapter  I 


INTRODUCTION 


The  problem  concerning  scattering  of  a plane  electromagnetic 
wave  from  a thin  conducting  square  plate  Is  easy  to  formulate  but  very 
difficult  to  solve.  The  major  difficulties  are  due  to  (1)  the  edge  con- 
dition, (2)  the  highly  if  not  nonintegrable  singularity  of  the  dyadic 
Green's  function  associated  with  a two-dimensional  surface,  and  (3)  the 
stability  problem  caused  by  its  thinness.  Obviously,  attempts  to  seek 
a closed  form  solution  are  bound  to  fail,  and  the  only  promising  approach 
Is  to  use  numerical  methods.  Among  those  who  did  work  In  this  area  using 
the  numerical  Integral  equation  solution  by  method  of  moments  (MOM)  are 
Ramat-Samii  and  Mittra[l],  and  Wang,  et  al  [2]. 


Although  the  moment  method  has  been  demonstrated  in  many  cases 
to  be  a very  successful  numerical  method  for  electromagnetic  problems, 
its  application  to  a two-dimensional  problem  has  not  been  good.  The  \40M 
offers  two  approaches  to  solve  the  surface  current  problems.  One  is  the 
wire  grid  modeling  which  suffers  from  many  numerical  and  physical  diffi- 
culties. The  geometrical  structure  Is  modeled  into  a wire  grid  repre- 
sentation and  the  results  are  sensitive  to  selected  wire  radius.  Grid 
representation  of  the  structure  has  loops  and  thus  loop  currents  can  be 
excited  which  are  manifestly  unphysical.  Again,  depending  on  the  grid 
structure,  internal  resonance  may  occur  and  thus,  alter  the  external 
resonances.  The  other  approach  is  the  surface  patch  modeling.  Here  the 
results  seem  to  vary  drastically  with  slight  changes  in  patch  size.  y,The 
matrices  tend  to  become  ill-conditioned.  Edges  present  particular  diffi- 
culties in  the  manipulation  of  the  integrand  in  the  patch  H modeling. 
Convergence  is  another  important  problem  in  this  approach.  In  order  to 
circumvent  these  drawbacks,  the  finite  element  method  [3]  is  used  in  this 
study.  The  finite  element  method  is  a relatively  old  approximation  tech- 
nique for  solving  boundary-value  problems  formulated  either  in  the  dif- 
ferential or  integral  form.  In  the  past,  this  method  had  been  mainly 
applied  to  the  field  of  structural  mechanics  [4]  and  only  recently  to 
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electromagnetic  problems  concerning  electrostatics  [5] [6]  and  propagation 
In  waveguides  [7] [8]. 

In  the  finite  element  method  approach  used  here,  the  solution  Is 
obtained  by  minimizing  a variational  functional  for  the  problem.  An 
approximation  to  a function  that  minimizes  the  variational  form  Is  con- 
structed from  combinations  of  certain  trial  functions.  These  trial  func- 
tions are  defined  on  the  region  In  which  a solution  to  the  differential 
or  Integral  equation  of  the  problem  is  sought.  The  region  is  divided  into 
subregions  which  are  called  elements.  Each  trial  function  is  zero  on  all 
parts  of  the  region  except  for  one  element.  The  various  trial  functions 
which  may  be  linear  functions,  polynomials  or'  the  like,  are  joined 
together  at  the  boundaries  of  the  element.  Values  of  the  trial  functions 
are  defined  at  certain  points  (nodes)  of  the  elements.  When  the  trial 
functions  which  are  also  called  shape  functions  are  substituted  in  the 
variational  form,  and  the  optimization  procedure  taken,  a matrix  equation 
will  result.  This  is  then  solved  to  determine  the  quantity  of  interest. 

i 

The  remaining  part  of  this  report  is  divided  into  five  chapters. 
Chapter  II  is  concerned  with  the  fundamental  concepts  of  the  finite  ele- 
l ment  method.  Chapter  III  presents  the  derivation  of  the  variational 

integral  equations  for  a conducting  flat-square  plate  and  a conducting 
bent-square  plate.  Chapter  IV  discusses  the  numerical  solution  of  the  var- 
iational integral  equations  by  the  finite  element  method.  Chapter  V 
presents  the  numerical  results  while  the  last  one.  Chapter  VI,  gives  the 
concluding  remarks  about  this  study.  Two  appendices  are  also  included  for 
completeness.  Appendix  A is  on  transformation  of  the  entire  surface  inte- 
gral into  a repeated  surface  integral  for  the  variational  integral  equation, 
and  Appendix  B gives  a brief  description  about  the  computer  programs. 

‘ 
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CHAPTER  II 

THE  FINITE  ELEMENT  METHOD  (FEM) 

2.1  Fundamentals  of  the  FEM 

The  FEM  is  primarily  a numerical  discretization  procedure  for 
solving  complex  boundary  value  problems.  The  method  was  originally  used  In 
the  field  of  structural  mechanics;  but  since  its  roots  belong  In  mathe- 
matics as  a class  of  approximation  procedure,  it  can  be  applied  to  a wide 
range  of  nonstructural  problems.  In  the  FEM,  the  region  of  the  problem 
is  divided  into  subdomains  or  finite  elements,  with  some  functional  repre- 
sentation of  the  solution  being  adopted  over  the  elements  so  that  the 
parameters  of  the  representation  become  unknowns  of  the  problem.  Usually 
the  element  parameters  are  the  nodal  values  and  their  derivatives  at  the 
nodes.  Although  the  region  of  the  problem  is  discretized  into  elements, 
the  whole  domain  remains  as  a continuum  because  of  the  imposed  restriction 
on  the  continuity  across  element  interfaces. 

The  FEM  itself  can  be  further  divided  according  to  the  procedure 
by  which  the  equations  in  the  nodal  values  are  formulated.  The  most  popu- 
lar method  is  the  so-called  variational  finite  element  method.  Here,  a 
stationary  principle  is  used  to  derive  the  set  of  equations  in  terms  of 
the  nodal  values  — often  a functional  obtained  from  the  variational  prin- 
ciple is  the  starting  point. 

2.2  Statement  of  the  Problem 

Let  L be  a linear  operator  defined  on  a dense  set  M of  the  com- 
plex  Hilbert  space  H with  inner  product  (u,v)  and  norm  (|u||  - (u,u)  where 
u,v£H.  In  general,  L has  to  be  neither  positive  definite  nor  symmetric 
and  it  could  be  either  differential  or  integral  or  integral-differential 
in  nature  [9].  If  the  governing  equation  is  represented  by 

Lu  - f , f C H (1) 

the  corresponding  stationary  functional  is  given  by 

F(u)  ■ (Lu,u)  - (u,f)  - (f,u).  (2) 
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Since  it  is  stationary  about  u,  if  the  approximation  solutions  differ  from 
the  exact  ones  by  an  amount  in  the  order  of  e where  e is  a small  quantity, 
then  the  calculation  of  F(u)  from  (2)  will  give  an  error  of  the  order  of 
e2.  Therefore  the  better  the  solutions  are  approximated,  the  closer  the 
results  will  be  to  the  true  values  of  u.  In  practice  the  stationary  prop- 
erty of  (2)  is  used  to  determine  the  approximate  solution  to  (1) . This  is 
called  classically  the  Rayleigh-Ritz  procedure  which  will  be  discussed  in 
the  following  sections. 

2.3  The  Subdivision  of  the  Region 

The  region  R is  subdivided  into  discrete  subregions  or  elements, 
each  of  the  same  general  form,  as  shown  in  Fig.  1.,  with  the  boundaries  of 
each  element  being  plane  or  curvilinear  faces,  and  with  the  adjacent  bound- 
aries of  any  pair  of  elements  being  coincident.  Commonly  used  elements  for 
surfaces  are  triangular  or  polygonal  form.  At  similar  positions  in  each 
element,  a number  of  points  are  identified  as  nodes.  They  are  generally  at 
the  vertices  of  the  elements,  and  at  positions  such  as  the  center  of  an 
edge,  the  centroid  of  a face  or  the  centroid  of  the  element  volume. 

Let  us  denote  the  nodal  values  of  the  solution  <(>  at  the  pc^  node  as 
Let  the  number  of  elements  into  which  region  R is  subdivided  be  Nc, 
and  the  total  number  of  nodes  in  R 1 D + B (Boundary)  be  n^  and  n^.  The 
total  number  of  nodes  in  a single  element  be  ns-  Then  the  nodal  values  of 
$ can  be  generally  expressed  as  a column  vector 


♦ i 

<f>2 


(3) 


k 


2-2 


r 


2.4  The  Element  Shape  Function 

To  solve  (2)  by  the  FEM,  one  needs  to  define  some  shape  functions 
or  Interpolation  functions.  These  functions  allow  us  to  express  the  solu- 
tion 4 at  any  position  of  R in  terms  of  only  the  nodal  values  {41}.  There- 
fore,we  assume  that  the  solution  $ an  be  prescribed  in  functional  forms, 
element  by  element,  across  the  region,  i.e.,  can  be  defined  piecewise  over 
the  region.  Within  each  element.  It  will  be  supposed  that  4 can  be 
described  by  a linear  combination  of  functions  Ni®,  N2®,  N. 


Ng®,  and  nodal  values  $j®,  $2®, 


k * 


$s  , thus 


<t>  =2Nie  + N2®  $2^  + n3C  $3*  + •••  N1, 


>ke  + •••  +Nse 


(4) 


or,  in  matrix  notation 
e „ e 


* =£  <Nie  N2e  ...  Nk®  ...  Nge)  {£*} 
e 

*Z(Ne) 

e 


(5) 

(6) 


Note  that  the  superscript  e is  used  here  to  identify  a particular  element. 

The  shape  functions  N®  are  restricted  to  being  functions  of  posi- 
tions. Since  the  true  solution  4 is  prescribed  as  being  continuous  and 
with  continuous  derivatives  (up  to  some  order)  across  the  region,  the 
piecewise  representation  (6)  should  have  the  same  properties.  Therefore 
the  shape  functions  are  restricted  by  the  following  conditions: 


1. 

2. 

3. 


N 


j 


■ 1 at  the  node 

■ 0 outside  element  e,  with  the  jth  node 

as  one  of  its  nodes 


■ a position  function  within  the  elements. 

In  choosing  the  shape  function,  one  has  to  pay  attention  to  convergence  in 
the  FEM.  Since  it  is  recognized  that  the  FEM  solution  to  a problem  with  a 
given  size  of  element  is  necessarily  an  approximation  to  the  exact  solu- 
tion, there  must  be  an  assurance  that  successive  finite  element  solutions 
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using  smaller  and  smaller  elements  will  converge  smoothly  to  the  exact 
solution  as  the  element  size  tends  to  zero.  While  comprehensive  condi- 
tions ensuring  convergence  are  not  yet  known  for  all  types  of  linear  prob- 
lems, there  are  certain  criteria  that  must  be  observed  in  order  to  obtain 
convergent  solutions: 

(1)  Completeness 

This  means  that  the  piecewise  representation  (6)  within  the  ele- 
ment of  the  variable/derivative  in  a key  Integral  must  be  capable  of  repre- 
senting any  continuous  function  as  the  element  size  decreases  to  zero. 
Mathematically,  the  piecewise  representation  calls  for  a complete  set  of 
functions  such  as  a polynomial  function  with  infinite  number  of  terms. 
However,  in  a FEM  representation,  only  a finite  number  of  terms  is  taken. 
But  as  pointed  out  by  Melosh  [10]  and  by  Zienkiewicz  [11],  a monotonic 
convergence  can  still  be  obtained  if  the  number  of  terms  used  in  the 
representation  allow  the  variable/derivative  up  to  and  including  order  t 
to  take  up  any  constant  value  within  the  element,  where  t being  the 
highest-order  derivative  of  the  variable  in  the  variational  functional. 

(2)  Compatibility 

This  means  that  the  representation  of  the  variable/derivative 
in  a key  Integral  of  (2)  must  tend  to  the  same  continuity  as  the  exact 
solution,  across  the  interelement  boundaries,  as  the  size  decreases  to 
zero.  If  for  a given  variational  functional,  the  highest-order  derivative 
involved  is  of  order  t,  the  derivatives  of  order  up  to  and  including 
(t  - 1)  are  known  as  the  principal  derivatives  of  that  variable.  Pre- 
suming that  the  exact  solutions  of  the  dependent  variables  are  continuous 
with  continuous  derivatives  up  to  at  least  order  t.  One  weak  requirement 
that  the  compatibility  criterion  is  satisfied  is  to  require  that  the  var- 
iable and  their  principal  derivatives  are  continuous  in  the  shape  func- 
tion representation.  This  means  that  the  highest-order  derivative  in  a 
key  Integral  will  have  a representation  that  is  at  worst  piecewise  con- 
tinuous, in  which  case  the  representation  will  tend  to  be  continuous  as  the 
element  size  tends  to  zero.  In  general,  completeness  and  compatibility  are 
sufficient  conditions  for  convergence  in  variational  finite  element  method. 


i 
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However,  these  conditions  are  very  strong  and  can  be  relaxed  [12],  In 
practice,  the  shape  functions  will  not  be  an  exact  representation  of  the 
true  solution,  but  an  approximate  one,  and  the  solution  obtained  will  be 
similarly  approximate. 

2.5  The  Subdivision  of  the  Functional 

Since  (1)  represents  essentially  a quadratic  function,  we  can 
write  it  as 

9 * f F(ui  ,U2,U3  ...  , ud  )dD  (7) 

JD 

where 

2 

F(ui  , U2  , U3  ...  , ud)  = an  u:  + a12  U1  u2  + al  3 U1  u3 

2 2 

+ a2i  u2  ux  + a.22  u2  +....+  a^d  ud  (8) 


D represents  the  domain  of  integration  which  can  be  a line,  surface  or 
volume,  and  ui,  u2»  U3  ...  , ud  represent  the  solution  $ and  its  various 
derivatives,  $x,  $y  ....  In  matrix  notation  (7)  becomes 


{u}T[A] 


{u}  dD 


where  [A]  is  a dxd  matrix  and  {u}  a-  dxl  vector,  or 


(9) 


[A] 


311  al2  * • ald 

3*21  a'22  • • • • a2(J 


adl  a-d2  ....  add 


{U} 


U1 

u2 


(10) 


(ID 
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and  superscript  T denotes  the  adjoint  of  a matrix.  In  general,  the  elements 
a^  are  functions  of  the  position. 

If  #e  is  the  contribution  of  an  element  to  the  total  integration 
in  (9),  then  (9)  can  be  written  as 


l 

E 

e“l 


* r 

&/ 


{u}1  [A]  {u}  dD 


(12) 


where  D represents  the  domain  of  element  e,  let  us  now  consider  a typical 
e 


term  uf  in  {u}  r ■ 0,  1,  2, 


By  definition,  uf  is  a spatial  deriv- 


ative of  $,  that  is  u ■ -2-  ♦ ■ D ♦,  where  S represents  a spatial  vari- 

r d§r 

able  of  concern. 


From  (6)  we  have 

u * (Ne)  {$e}  in  element  e. 

Thus  within  element  e 

ur  - Dr4  - (Dr  Ne)(te}  - (Ure){*e) 

when  (u  *)  represents  the  row  vector  for  the  r-derivative  of  the  shape 
function.  So  applying  (13)  for  every  element  we  obtain 

{u>  - (U)  {**} 

where 


(13) 


(14) 


[U]  - 


Di  N!6 

Di  N2 

D2  NX6 

• 

d2  n2 
• 

• 

• 

Dd  NXe 

• 

• 

Dd  N2 

D2  N 


D . N„ 
d s 


(15) 


Substitution  of  (14)  into  (9)  yields 


* “ 2Z  / {*e>T[u]T[A][0]Ue>  dDe 


e - I D. 


which  shows  that  $ is  now  a function  of  the  n^  nodal  values  $2»  •••  • 


2.6  The  Stationary  Condition 

In  order  to  solve  (16)  we  have  to  invoke  the  variational  prin- 
ciple. The  condition  that  * is  stationary  is  given  by 


hi  _ hi-  m „ 


3 i/d  <Jij 

8 $/3  ^2 


3 i/d  i 


From  (12)  we  get 

t e 

£ tfcr0 

e ■ 1 

2 . 7 The  Element  Matrix  Equation 


To  get  the  element  matrix  equation  we  have  to  combine  (19)  and 
(16) . Considering  the  term  f°r  an  element  e » we  ®et 


- = f —2—  [ue}T[U]T[A]  [U]  {<pe}l  dDe 
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Ml_  - Ke]^e> 

d(*e} 


(25) 


where 


[A1*]-/  2[Al6]dDc 

^ T) 


(26) 


If  we  substitute  (15)  into  (20)  and  carry  out  the  mathematics  we  will 
■ th  •> *-  r »e- 


obtain  for  the  ijth  element  of  [B  } as 


'ij 


' J ' 1 |^DlNie(an  °1  Nj6  + a12  D2  Nj6  +...•  + aid  Dd  N^e) 


+ Di:Nie(a2i  Dj  N e + a22  D2  N 6 + 


+ a2d  Dd  Nj  ) 


+ Dj  N1e(ad]  Dj  Nj6  + ad2  D2  Nj  + 


+ a 


dd 


»„  »/>] 


dD  (27) 
e 


Note  that  in  (27)  the  subscripts  on  the  Ne  are  in  terms  of  the  node  ident- 
ifiers, not  system  node  numbers. 

Note  that  the  shape  functions  are  explicitly  defined  functions 
of  spatial  variables.  The  Integrand  of  a particular  term,  say 


/ d0. 


could  be  evaluated  as  an  explicit  function  of  x,  y and  z.  If  a^  are 


constant  coefficients,  the  prescribed  Integration  over  the  defined  domain 
of  the  element  would,  in  consequence,  evaluate  the  term  as  a scalar. 
The  Integration,  if  simple,  could  be  carried  out  analytically.  However, 


if  a^j  are  complex  functions  of  x,  y and  z,  then  the  integration  would 


generally  require  a numerical  solution.  Therefore,  the  computational 
time  involved  in  £ 
complex  functions. 


time  involved  in  a problem  depends  very  much  on  whether  a^  are  simple  or 
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r 
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2.8  The  Boundary  Condition  (B.C.) 

It  Is  known  in  boundary-value  problems  that  the  solution  is  not 
unique  unless  it  meets  all  the  required  boundary  conditions.  However,  in 
the  variational  finite  element  methods,  if  the  specified  boundary  condi- 
tions are  natural  boundary  conditions  for  the  problem,  then  it  can  be 
shown  that  the  class  of  admissible  functions  is  not  required  to  satisfy 
these.  In  order  to  illustrate  the  treatment  of  the  boundary  condition  in 
the  matrix  equation  (25)  let  us  assume  a Dirichlet  boundary  condition  such 
that 

$ - g(x,y,z)  on  B.  (28) 

Using  (28),  the  n^  nodal  values  ($p)jj  for  the  boundary  nodes  on  B can  be 

calculated  yielding  n.  equations  of  the  form 


(29) 


« 


* 


which  implies  that  if 
a constant  value,  then 


satisfies  the  boundary  condition  and  hence  it  is 

■ 0 for  an  element  containing  node  p. 

P 


Thus  to  include  the  B.C.  in  the  element  matrix  equation,  the  simplest  pro- 
cedure is  to  replace  the  pc^  row  of  the  matrix  [A1®]  in  (25)  by  the  row 
matrix  of  (29).  In  other  words,  if  p is  a boundary-condition  node,  put 
zeros  in  the  pc^  row  of  the  fAie]  in  (25)  except  for  a 1 in  the  diagonal 
position  and  put  in  the  pc^  row  of  the  driving  vector  the  boundary  value 
given  by  (28). 
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CHAPTER  III 


DERIVATION  OF  VARIATIONAL  INTEGRAL  EQUATIONS  FOR 
A CONDUCTING  FLAT-SQUARE  PLATE  AND  A CONDUCTING 
BENT-SQUARE  PLATE 


3.1  Flat-Square  Plate 

We  are  considering  a thin  (compared  to  the  wavelength  of  the 
Incident  wave),  perfectly  conducting  flat-square  plate  as  shown  In  Fig.  2. 
The  plate  Is  Illuminated  by  a plane  electromagnetic  wave  with  Its  E vec- 
tor polarized  at  an  angle  8 with  respect  to  the  positive  x-axls,  and  Its 
propagation  vector  k lying  on  the  y-z  plane  at  an  angle  $ with  respect  to 
the  negative  y-axis.  The  thickness  of  the  plate  Is  2t  and  the  dimensions 
of  the  plate  are  2a  x 2b.  The  MRS  units  and  e^U)t  time  variation  are 
adopted.  It  is  the  objective  of  this  study  to  compute  the  induced  surface 
current  density  due  to  an  Incident  plane  wave  by  the  FEM. 


in 


Let 


E(r)  - (cos  6 x + sin  6 y)  co*  ♦ ’ ' 8ln  ♦> 


(30) 


The  E-field  integral  equation  governing  the  Induced  surface  current  Is 
in 


E(r)  “ - jwp0  r J(r»)  • G(r|r')dS' 
''s' 


(31) 


where  £,r'  are  the  observation  and  source  points  on  the  plate  surface 
respectively, 

in 

E(r)  ■ the  incident  electric  field  (v/  ) 

— — m 

J(r')  ■ the  surface  current  density  (a/m2) 

u ■ angular  frequency  (rad/sec) 

W0  ■ the  permeability  of  free  space  (**/  ) 

j - /=T 

S’  ■ the  surface  area  of  plate  (m2) 

G(£|r*)  ■ the  free  space  dyadic  Green's  function 

- (i  + ^2  sClIi') 


(32) 
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Figure  2.  Geometry  of  the  Flat-Square  Plate 
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0 


I ■ the  unit  dyad  and 


i(ik')  - «o  " 


[!«' 


2 + t2 


Kouyoumjian  [13]  has  shown  that  the  following  is  a stationary  quantity 


■ -2  f E^U(r)  • J(r)dS  + juyo  f f J(r)  • G(r|r ')  • J(r')dS'< 

C *C  * /o 


For  a two-dimensional  problem  G(_r  | ’ ) takes  the  following  form 


G(r|r')  - 


i a?  - 

i + —j  X 

k2  dx2 


1_  d 

V2  dxdy 


XxxX  Xxy  y 


X x X y 
yx  yy  , 


i_  a2  - 

k2  axay  y 


i + 1 — * 

T2  a^y 

k °y 


g(r|r') 


where 


5s  [(«• + i)  - 0/‘]  - 0 


X - xw 
xy  yx 


U-x’)(y-y')  [(q  + r)  -r] 


-r  V1’ [(“♦*)* -*]' 8 


L 


with 


8 - Jk  + r 


(40) 


Letting 

J(r)  - Jx(£)x  + Jy (r)y 

and  substituting  (36)  and  (41)  into  (34)  we  obtain 


(41) 


- -2  f E*  (r)  • J(r)dS  + ju>P0  f f \ ^^(r)  Jx(r  ’ ) 

+ Vi>V£'>  + Xxy  > + Xyy  >]dS'dS 


(42) 


(42)  represents  the  variational  integral  equation  for  the  flat-square  plate. 
3.2  Bent-Square  Plate 

Here  we  are  concerned  with  a bent-square  plate  as  shown  in  Fig.  3. 
The  plate  is  bent  at  an  angle  <^.  Because  the  problem  is  symmetrical 
about  the  y axis  we  can  set  up  the  x-y-z  coordinate  system  with  its  origin 
on  the  y-axis  and  one  of  the  edges.  For  mathematical  convenience,  we  also 
set  up  another  coordinate  system  x^-y^-^  on  t*ie  Inclined  surface  and  the 
bent. 

Now  for  a three-dimensional  problem,  the  free  space  Green's 
function  takes  the  following  form 


G(r|r’) 


1 

_di 

A A 

XX 

i a2 

i_  a2 

k2 

ax2 

k2  axdy  ^ 

k2  axdy 

1 

a2 

A A 

, . i a2  ~ 

i a2 

k2 

aya» 

yx 

i+? 

k2  ayaz 

1 

i 

A A 

ZX 

i a2  ** 

— ■'  ■!  — ■ 7Y 

i + JL 
k2  az2 

k2 

aza» 

k2  azay 

yz 


zz 


g(r|r’)  (4: 
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r 


XX 

X 

xy 

Xv„  xz 

XX 

xy 

xz 

A A 

A 

A A 

A A A 

yx 

yx 

A 

yy 

yy 

> 

N 

N 

A A 

A A 

\ A A 

zx 

ZX 

X 

zy 

zy 

Azz 

L ' S 

are 

similar  to 

those 

given  : 

(44) 


We  now  divide  the  surface  area  of  the  plate  into  two  regions:  one 


is  called  the  bent  or  inclined  region  denoted  by  and  the  other  one  is 


called  the  unbent  region  denoted  by  Su.  Because  of  this  we  can  divide 


the  area  of  integration  into  four  parts: 


(1)  when  r and  r'  are  both  in  S : 
— — u 


Here  we  let 


J(r)  = J (r)x  + J (r)y 

x y 


(45) 


J(r')  - Jx(r»)£  + Jy(r’)y 


Substituting  of  (45)  into  (34)  and  using  (44)  we  get 


—2 'f  E(r)  • J(r)dS  + ju)p0  f f {"^^(r) Jx(r ' 

\ A.  J s: 


u 


u u 


+ V^)Jy(l,)  + XyxJy(-)Jx(-'}  + xyyJy(-)Jy(~  }]  ds'  ds  (46) 


where 


(r,r')CS 


(2)  when  £ and  r/  are  both  in 
Here  we  let 
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J 


r 


i 


J(r)  - J^Cr)^  + Jyb(r)yb  - Jx>^(r)i  + cos  p Jyh(r)y  + sin  p Jy^(r)z 


*b 


yb 


yb 


(47) 


j(r')  “ Jxb(E.' )*b  + Jybfe'^b  " Jxb^-'^x+  cos  p Jyb^’)y  + 8ln  p Jyb^-' ^ 


where  p ■ 180°  - $b>  subscript  b refers  to  the  x^)_yb-zb  coordinates  and 
from  (34)  we  get 


*2  " -if  Exn(r)  JXfe(r)dS  + jmp0 ^ y 


cos  p Jyb  y + sin  p J. 


W 


yb(£)J] 


jJ*b(- + V + V]  + coa  “ Jyb<£'  > [xy**  + v + v] 


+ 8l»  » Vi'1  [Az,;  + xb/  + Xj]  j dS'dS 


(48) 


Since  x - 


y - cosp  yfc  - sin  p zfc 


(49) 


z - sin  p y,  + cos  pz. 


b - '•w“  K‘b 
therefore  (48)  reduces  to 


*2 


‘-if  E(?)  • J(r)  dS  + jtoy0  f f ( X J (r)  J (r* 

Vvl 


b b 


[C°8  P Ayx  + 8ln  P XZX]  \ <^>  V£,)  + (Xjjy  COS  p + Xxz  sin  p) 

yb(l)  JXb(x’)  + [COS  P (Xyy  COS  P + Xyz  8±n  P)  + 8in  P (XZy  008  P 

P)]  Jyb(£)  Jyb(r')|  dS'dS. 


+ Xzz  8ln 


(50) 


! 
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Also  since  x 


y - yb  cos  p + aj 
* ■ y^  sin  p 


and  R - [(x-x')2  + (y-y')2  + (z-z')2  + t2] 1/2 

■ [(xb_xb,>  +(yb  cos  p-yb'  cos  o)  +(yb  sin  p~yb'  sin  p) 2+t 2J 
“ [(xb  " V)2  + <yb  “ yb,)2  + t2]V2  (52) 


dR  _ dR 

d**,  " d* 

d R _ dR  dy  , dR  dz 
^yb  dy  ^yb  *z  dyb 

- cos  0 + sin  „ |S 


(co* “ h + sl"  p f^)  (“8  s ^ 


+ sin  p 


cos2  P ^-^5  + sin  p cos  p 

d y2  dzd’ 


^2  2 2 

■r~—  + sin  p cos  p ^ — ■ + sin2  p 
dzdy  dyd*  dz2 


Combining  (53)  to  (55)  we  have 


cos  p A^  + sln  p Xzx  . ^ 
cos  s Xxy  + sln  » x*z  ' V» 


cos  p (*  cos  p + A sin  p)  + sin  p (A  cos  p + A sin  p)  - A 


ybyb  (56) 


and  so  (50)  becomes 


J(r)dS  + J<uu0 


*ss  S® 


J_  (r') 
*b 


+ X, 


ybxb 


s® 


Jybl- 


) + V 


Vb 


V£,) 


Jyb<r) 


+ X 


vb 


V 


where  (^r'JCS^ 

(3)  when  £ is  in  Su  and  r'  is  in 
Here  we  have 

J (r)  -Jx  (r)x  + Jy  (r)y 

J(r')  - + Jyb^-^b  ~ Jxb(l.'^  + co8  p Jyb(l')y  + 

Substitution  of  (58)  into  (34)  yields 


r in 
-2  / E(r) 

*/s 


J(r)  dS  + ju»u 


+ (cos  p Xyx  + sin  p Xzx)  Jx(r)  Jyb(r')  + (cos  P ^ + sin  p Xzy 
* \<S')  | ISMS 

where  r CSu;  r'csfe 


(4)  when  r is  in  S.  and  r'  is  in  S 

D — U 

Here  we  have 

J(r)  - JXb  (r)  ^ + Jyb(r)yb  - JXb^)*  + cos  P Jyb^)y  + sin  p Jy^  (r ) z 
J(r')  - Jx(r')x  + Jy(r')y 
Substitution  of  (60)  in  (34)  gives 


(57) 

r')jdS’dS 


(58) 

sin  p Jyb(£’  )* 


Jy(r)  Jyb(*f) 
(59) 


(60) 


1(1)  • J(l)<*S  + Juvo J J j^Jj^Cl)  Jx(r’) 

+ *yx  ^(1)  Jy(l')  + (cos  p >xy  + sin  p *xz)  Jx(l’)  Jyb(r) 

(61) 

+ (Xyy  cos  p Xy^  sin  p)  Jyb(l)  Jy(r' )|  ds'ds 

where  r £S.  ; r'  £ S 
— b—u 

Equations  (46),  (57),  (59)  and  (61)  constitute  the  variational  integral  equa- 
tions for  the  bent-square  plate.  These  will  be  solved  by  the  FEM  in  the  next 
chapter. 

3.3  The  Thickness  Parameter 

In  our  solution  of  the  problem,  we  have  introduced  a thickness 
parameter.  The  primary  reason  for  doing  this  is  to  circumvent  the  difficulty 
of  treating  the  singularity  of  the  E-field  integral  equation.  This  singular- 
ity, if  without  any  modification,  is  considered  not  integrable  for  the  two- 
dimensional  case,  although  there  exists  some  analytical  approach  for  the 
three-dimensional  case  [14] . Physically,  this  introduction  of  a thickness 
parameter  is  equivalent  to  using  the  so-called  extended  boundary  condition 
method  [15]  in  which  it  is  required  that  the  incident  field  be  cancelled 
by  the  scattered  field  in  an  arbitrary  region  completely  within  the  scat- 
terer,  so  that  the  field  point  and  the  observation  point  would  never  coin- 
cide. One  disadvantage  of  employing  a thickness  parameter  is  that  instead 
of  considering  a two-dimensional  problem,  we  are  actually  solving  a three- 
dimensional  problem.  In  order  to  restrict  ourselves  to  a two-dimensional 
problem,  we  have  to  assume  that  the  thickness  of  the  plate  has  to  be  small 
compared  with  the  wavelength  of  the  incident  field  and  the  dimensions  of 
the  plate. 
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CHAflfiR  iV 

NUMERICAL  SOLUTION  OF  THE  VARIATIONAL  INTEGRAL 
EQUATIONS  BY  THE  FINITE  ELEMENT  METHOD 


In  this  chapter  we  discuss  the  numerical  procedure  of  solving 
the  variational  Integral  equations  developed  in  Chapter  3 for  both  the 
flat-  and  bent-square  plates  using  the  FEM. 

4.1  Discretization  of  the  Surface 


As  shown  in  Fig.  4 the  square  surface  is  subdivided  into  triangu- 
lar elements.  In  the  discretization,  consideration  is  given  to  the  struc- 
tural symmetry  with  respect  to  the  incident  field.  For  the  case  of  normal 
incidence,  there  are  two-fold  symmetry  about  the  x and  y axes,  hence  it  is 
sufficient  to  represent  the  unknown  nodal  values  over  only  one  quadrant  of 
the  surface.  For  oblique  incidence,  we  have  just  one-fold  symmetry  about 
the  y-axis;  therefore,  we  have  to  represent  the  unknown  nodal  values  over 
one-half  of  the  surface.  The  contributions  due  to  the  other  parts  of  the 
surface  in  computing  the  matrix  elements  can  be  taken  care  of  using  the 
image  techniques.  Once  the  surface  has  been  subdivided,  system  node  num- 
bers and  element  numbers  are  assigned  to  each  element.  For  example,  for 
the  normal  incidence  case,  as  illustrated  in  Fig.  4 (a),  the  system  node 

numbers  are  from  1 to  16  ( Nn  * n , + n.  * 16 

/ t d \ b 

ments  for  one  quadrant  is  18  IN  * 18). 


^ and  the  total  number  of  ele- 


Shape  Functions  and  Area  Coordinates 


In  order  to  compute  the  matrix  elements,  it  is  necessary  to  carry 
out  the  integration  over  the  elements.  To  do  this  we  have  to  generate  some 
interpolation  functions  or  shape  functions.  The  generation  of  interpola- 
tion functions  for  a triangular  element  is  simplified  considerably  if  one 
works  with  area  coordinates.  As  shown  in  Fig.  5,  the  vertices  of  the 
triangle  are  numbered  by  £,  m and  n.  The  partial  area  opposite  to  node  t 
is  denoted  by  A^,  and  its  coordinates  in  the  x-y  coordinate  system  are 

(*£,  The  total  area  of  the  element  is  designated  by  A.  For  any  point 

in  the  triangle  we  assign  three  coordinates  denoted  by  L,,  L and  L . 

•cm  n 
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Subdivision  of  Che  Surface  inco  Triangular 
Eleaencs.  (a)  for  Normal  Incidence  (b) 
for  Oblique  Incidence. 


(67) 


Note  that  only  two  of  the  area  coordinates  are  independent  because 

L»  + L + L - 1 
t m n 

Since  the  expressions  for  the  matrix  elements  involve  Cartesian  deriva- 
tives, it  will  be  necessary  to  transform  the  derivatives  too.  The  trans- 
formation laws  are  readily  deduced  from  (65) 


3x  [f(Li  ’ Lm  » Ln)] 


T 3JL3Jl 

iTl  3Li  3x 


f y li- 
fe i3Li 


(68) 


37  tf(Li  ’ Lm  ’ Ln» 


y x — 
ft  i3Li 


Higher  derivatives  are  generated  by  repeated  application  of  (68) . For 
example. 


2 3 3 2 

3 f £ £ ^ — 


3x3 y ft  ft.  1 j 3Li3Lj 


(69) 


For  Integration  ve  can  easily  show  that 


/ 


Vs ' t 


(70) 


/L.L  dS  = 7T  2 1*  m 
2 m 12 


(71) 


/ 


<I.t,2dS  - f 


(72) 
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and  the  following  general  formula: 


( 1 ! ) ( j ! ) ( k ! ) 
(i+j+k+2) ! 


2A 


where  i,  j and  k are  positive  integers. 


(73) 


In  the  above  analysis,  we  have  assumed  a linear  interpolation  function  or 
a first-order  polynomial  for  the  shape  function  since  there  are  totally 
only  three  nodes  available  in  an  element.  However,  for  higher  accuracy, 
better  approximation  or  for  convergence  reason,  higher  order  polynomials 
may  be  used.  This  can  be  done  by  using  more  nodes  by  placing  them  on  the 
sides  of  the  triangles  or  inside  the  triangles  as  interior  nodes.  Another 
shape  function  that  can  easily  be  generated  from  the  area  coordinates  is 
th4  sinusoidal  shape  function  which  is  defined  as 


sin 


(¥) 


and  L 


(74) 


4.3  Conversion  of  the  Integral  Equations  into  Matrix  Equations 

Let  us  focus  our  attention  on  the  integral  equation  for  the  flat- 
square  plate.  The  integral  equation  for  this  case  is  (42). 


s s' 


+ X J (r) J (r')  + X J (r ' ) J (r)  + X J (r')J  (r)]dS'dS 
xy  x — y — xy  x — y — yy  y — y — 


In  terms  of  the  discrete  triangular  elements,  and  using  i and  j to  ident- 
ify the  unprimed  and  primed  coordinates  respectively,  the  above  equation 
becomes 


where  A^  and  Aj  represent  the  triangular  elements  for  the  observation 
point  and  the  source  point  respectively.  Now  we  can  express  the  current 
density  J_  in  an  element  in  terms  of  their  nodal  values.  Thus  for  the  x 
components 


li  mi  ni 

J (r)  • L„  J + L J + L J 
x — x x ni  x 


“i  "j 

J (r')  =>  L„  J J + L J 3 + L J J 

x — x nij  x nj  X 

and  similarly  for  the  y components 

*i  mi  ni 

J (r)  * L.  J + L J + L J 

y - y m y ni  y 


S mi 

J (r')  = L.  J ^ + L Jj+L  Jj 

y - y “j  y nj  y 


(76) 


(77) 


Now  let  us  write  (30)  as 


in 


E(r)  - E (r)  x + E (r ) y 


in  . ^ 


where 


in  -jk(y  cos  $ - z sin  <fr) 

Ex(r)  - cos  0 e 

in  -jk(y  cos  $ - z sin  $) 

E (r)  ■ sin  0 e 

y- 


Substitution  of  (76)  and  (77)  into  (75)  gives 


4-7 


r 
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* - -2  e / \c{Lhj!1 + v-"1 + v*"1)*  ^(v*  1 + v»  1 

4i 

+ v/j)  as  ♦ i»u0  i;  % f f jv^  + + V-"1) 


(v* 3 + V* 3 + \3* 3) 

+4t/:^sC^sO)(v/l+vy“i+sjy"1) 

* + V ' + + SJyl)  + V1) 

(v/1 + v/1 + V'Kv1 + V*”1 + Vy”1ldS' 


(78) 


+ A 

yy 

where 


dS 


1^  , ■ 1 • 2 N 


i , 3 » 1 . 2 


end 


k k 

N11  “ total  number  of  nodes  for  one  quadrant.  J and  J are 

. * y 

the  unknown  nodal  values  for  J and  J at  the  k^“  node.  To  solve  for 

if  ir  ^ y 

J * and  J K we  differentiate  (78)  with  respect  to  J and  J and  set 
x y * y 

equal  to  zero  since  • Is  stationary;  thus,  we  get 
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-^T"2  L [ *tn(h  6 1 + Lm  6m  + Ln  6n  W 

3J  k i?l  J X \ *i  *ik  “i  “ik  "i  "Ik/ 

Ai 

4NC  4Nt  C fi  \i  \ 

+ L E / / )X  L 0 6 0 +L  6 +L  « 1 

° i-i  j-i  y y j xx|\  *i  fcik  mi  mik  ni  nik/ 

(L  J*  ^ + L J J+L  J /l„  <5„  + L 6 + L <5  \ 

li  X °j  X "j  X / \ V mj  mjk  nj  njkj 

/l.  Ji  + LJ1  + LJ1\+x/L„6„  +L6  +L5  \ 

\ *i  X “i  X ni  X !\  Xy\  *3  3k  “j  mjk  "j  njkJ 

(l.  J * + L J A + L J i)+  X (l.  6.  + L 6 + L 6 \ 

\*iy  mi  y niy  / xy\£iAik  mi  mik  ninikj 

(h'y*  + + L„jJy”J)|^'dS  - 0 


where 


6.  - 1 if  l . - k 

*ik  1 

- 0 otherwise. 


In  a similar  manner  we  can  obtain  — — — which  is  deducible  from  (79)  by 

3J  k 

y 

replacing  X by  X , E*n  by  Ein  and  J by  J . 

xx  1 yy  x ; y x 1 y 

In  matrix  notation  we  have 

[S]{J>  - {E>  (8i 

where  [Sj  is  a ( 2Nn  x 2Nn)  "Impedance"  matrix 

{J}  is  a (2Nn  x l)  column  vector  for  the  unknown  nodal 

current  densities. 

(E)  is  a (2Nn  x l)  column  vector  for  the  driving  incident 

field. 


L 


and  the  matrix  elements  are  computed  from 


4Nt  4NC 


S,  ■ imp 
kq  J o 


4N  4N  r r ,,  , 

E E //  ^ + V 

i-1  j-1  JJ  **  mn  1 *J  *jq  mj 

A A ' 


Ai  AJ 


+ /l.  6,  + L 6 + L 6 Yl  dS’ 

d l,  l.  m.  m.  n.  n.  II 

mn\  i iq  1 lq  1 lql 


dS 


for  1 < k , q < N 


4NC  4NC 


Sk,  ‘ *"*0 


WH  fMJb 


X 16  *"IL.  6.  + L 6 

xy  d l d.  d.  m.  m 

1 ran\  j jq  j 


4i  4i 


+ [l,  8,  + L 6 + L 6 \ 

d I d.  d.  m.  m.  n n.  I 
mn\  i iq  1 iq  i iq  / 


dS'dS 


for 


1 ^ k < Nn  and  N° . _<  q + Nn  < M 


Nn  <_  k + Nn  £ M and  1 < <i  < N™ 


rt 

4Nt 

Nn  " 

jwv  £ 

E 

i-1 

i=l 

6jk  | 

/ 

L„  6.  + 

L 

4mn' 

d.  d. 

V i iq 

mi 

A1  *j 


jq 


'dS 


for 


Nn  _<  k + Nn  , q + Nn  M 


+ L 6 
m.  m. 


+ L 


6 

n,u 


where 


9 


As  an  example,  let  us  consider  Fig.  4(a).  We  have  here  a total  of  16 
nodes.  Therefore,  we  should  have  16  unknown  nodal  values  for  Jx  and 
another  16  for  Jy.  Hence,  the  matrix  size  will  be  32  x 32.  The  first  16 
rows  of  the  matrix  represent  Jx  coupling  to  Jx  and  Jy  whereas  the  next  16 
rows  are  for  Jy  coupling  to  Jx  and  Jy.  For  a typical  node,  say  node  2, 
we  get  from  (84)  to  (86) 


S - 4 jaip0  \ / f + f f + f f \ L L dS'  dS 

* \ k \\  W 2 2 

The  factor  4 accounts  for  four  quadrants 

- ■ >-\( H (A kl H l 

A1  A1  A1  A2  A1  A7  A2  A1  A2  A2 


+//+/  / +/  f *f  f | 


UL  L dS'  dS 
2 5 

+ Images 

Contributions 


18  18 


■ 4 jam- 


WJAA,  I 


Xyy  ^ ^ dS'  dS 


and  so  on. 


4.4  Numerical  Integration 

By  using  the  FEM,  the  integration  over  the  whole  surface  is  now 
reduced  to  a summation  of  integrations  over  the  individual  elements.  The 
integration  over  each  element  can  be  carried  out  numerically  using  one  of 
the  existing  numerical  quadrature  formulas  such  as  Simpson's  rule,  Gauss- 
ian quadrature,  Romberg's  extrapolation  method,  etc.  However,  the  simplest 
method  is  to  replace  the  integration  by  its  Riemann  sum  with  unit-weighting 
coefficients.  That  is,  if  we  divide  the  ic^  element  into  N subareas  we 


F 


where 


/« 

f(x,y)  dxdy  f^xij*yij^  ASij 

«i  j - 1 

■ the  area  of  the  element 

ASjj  ■ the  j1*1  subarea  of  the  ic^  element 

xij  *ylj  " t*ie  X an<*  y coor^inates  °f  the  centers  of 
the  subarea  of  the  i**1  element. 


Integration  over  the  Singular  Element 

Since  the  singular  element  involves  singularity,  the  Integra- 
tion over  this  element  should  be  performed  more  accurately  than  the  non- 
singular (regular)  element.  As  shown  In  Fig.  6(a),  the  singular  element 
is  divided  equally  into  square  and  isosceles  triangle  subdivisions  (size 
6fl).  Within  each  subdivision,  the  current  is  assumed  constant  and  repre- 
sented by  the  value  at  the  center  of  the  square  or  the  midpoint  of  the 
diagonal  of  the  triangle.  If  the  field  point  and  the  source  point  are  in 
the  same  subdivision,  then  the  subdivision  is  again  divided  into  smaller 
cells  (size  6SS).  In  this  way  the  accuracy  of  integration  over  the  sin- 
gular element  can  be  controlled  by  specifying  two  numbers  — one  for  the 
division  of  the  singular  element  and  the  other  for  the  subdivision.  It 
should  be  noted  that  the  integration  over  the  singular  element  is  per- 
formed only  once  in  the  program. 


Integration  over  the  Regular  Element 

Since  less  accuracy  is  usually  required  for  integration  over  a 
regular  element,  the  regular  element  is  normally  divided  into  a larger 
subdivision  (size  fi)  as  shown  in  Fig.  6(b).  No  further  dividing  is 
needed  for  the  subdivision. 
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4.5  Symmetry  Consideration 

Although  structural  symmetry  allows  us  to  reduce  the  total  nodal 
unknowns  of  the  problem  to  those  for  either  one  quadrant  or  one-half  of 
the  surface,  in  computing  the  matrix  elements  we  still  have  to  carry  out 
the  integration  over  the  whole  surface.  However,  by  using  the  image 
technique,  we  can  restrict  our  integration  to  only  either  one  quadrant  or 
one-half  of  the  surface.  This  is  feasible  because  of  the  following 
consideration: 

Without  loss  of  generality  let  us  use  (84)  to  calculate  S^.  It  is 
seen  that  every  term  in  the  integrand  is  constant  except  X^  which  is  a 
function  of  positions  for  the  observation  point  and  the  source  point. 
Moreover,  for  each  element  in  the  first  quadrant,  there  will  be  a similar 
element  (image)  in  the  other  quadrants  or  other  half.  Therefore,  assuming 
normal  incidence,  we  have  three  images  (Fig.  7)  and  if  we  place  the  obser- 
vation point  in  the  first  quadrant  we  can  write  (85)  as 


Skq  = j4“Po 


am 


2 3 4 

+ X + X + X 

XX  XX  XX  XX 


)K( 


L*  x 

j jq 


+ L 6 + L 5n  ^ + 5lk  SA.  + V,  V,  + Lq4  &nia)  dS 

“j  “jq  nj  njq/  M ^ i iq  i 


(92) 

dS 


4 ■ '<*’  • y') 

(first  quadrant) 

4 ■ s<*+  • 50 

(second  quadrant) 

^ . + +. 
X = §(x  , y ) 

XX 

(third  quadrant) 

4 _ + 

X * §(x  , y ) 

XX 

(fourth  quadrant) 

§<u  • v) 

(93) 


(94) 
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4.6  Integration  over  One  Side  of  the  Plate 

In  the  surface  Integration  of  (34),  the  Integration  has  to  be 
carried  out  over  the  entire  surface  of  the  plate.  If  vre  neglect  the  con- 
tribution due  to  the  sides  corresponding  to  the  thickness  (this  is  not  a 
bad  approximation  since  the  thickness  is  small  compared  to  the  other  sides 
and  the  wavelength  of  the  incident  field)  and  represent  surface  areas  of 
the  top  side  and  the  bottom  side  by  S+  and  S_,  we  can  write  (34)  as 


* - -2 


Si 


ln<I>  * J(r)dS  - 2 j E±n(r)  • J(r)dS 


+ juu  ) f [ J(r)  * G(x'x')  • J(r')dS' 

k s; 


dS 


(97) 


+ II  J(r)  • G(r’r’)  • JCr^dS'dS  + / / J(r)  • G^'r')  • J(r')dS’ 

+ //  J(r)  • G(r'r')  • J(r’)dS'dSj 
s_  S_  J 


dS 
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It  is  known  that  for  sufficiently  large  plate  sizes  the  current  on  the  bot- 
tom side  S_,  which  is  not  facing  the  incident  field,  is  relatively  small 
compared  to  that  of  the  top  side  S+  except  near  the  edges  where  the  current 
exhibits  singular  behavior  (for  straight  edges,  the  current  component 
parallel  to  the  edge  behaves  like  while  the  component  perpendicular  to 
the  edge  behaves  like  /s  » where  s is  a small  distance  from  the  edge) . 
Therefore  we  can  neglect  the  integration  over  S except  the  small  region 
close  to  the  edges.  Let  us  denote  that  small  region  by  6S  , and  so  (94) 
reduces  to 


■V 

^n(r)  • J(r)dS  + jtop 

o 

j // 

s+  +<5S_ 

(s+  s; 

• G(r'r')  • J(r')dS'dS 

*// 

s+  6S1 

6s_  s'+ 

\ 

*f  f±(i) 

<5S_  6S'_ 

• G(r’£')  • J((_r')dS'ds| 

J(r)  • 2(r'r')  • J(r')dS'dS 


i’dS  (98) 


which  can  be  combined  into  the  following  form 

* m ~2f  y J(r)  • G(r ' r' ) • J(r' 

§+  +<5S_  / s+  -H5S_  s| 

f J(r)  • G(r'jr)  • J (j:' 


)dS'dS 


)dS'dS) 


(99) 


S+  +<5S  6S' 


or 


/•in  r _ 

* - -2  / E(r)  • J(r)  dS  + Jwp0  / f J(r)  • G(r|r')  • J 

*S+  + «S_  S+  + «S_  *4+’  + 6S^ 


J(r)  dS'  dS 
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If  we  let  d(s  - 6S_)  denote  the  current  density  on  6S_,  where 


«(S  - «S_J  - a value  if  r C6S 


■ 0 otherwise 

and  note  that  the  current  on  6S_  has  a reverse  direction,  then  we  can 
write  (96)  as 


♦ 


J(r)  [ 1 - 6(8  - «Sj]  dS 


(100) 


+ J“U0 


s s 


«(s  - 6S_)J  J(r) 


G(r|r')  • J(r')[  1 - «($'  - «S^)]  dS  dS 

(101) 


Equation  (101)  now  shows  that  the  integration  is  restricted  only  to  the  top 
side  of  the  plate.  Note  that  the  factor  6(s  - 6S_)  is  equal  to  zero  when 
the  observation  point  is  far  from  the  edge.  When  the  field  point  is  close 
to  the  edge,  d(S  - dS  ) can  take  on  any  values  larger  or  smaller  than 
zero.  From  study  of  canonical  problems  concerning  currents  near  edges  of 
thin  objects,  we  can  approximate  6(s  - 6S  ) by  the  following: 


For  the  parallel  component 


«(s  - 6S_)  - 1/2 

when  £ r 6S_ 

(102) 

For  the  perpendicular  component 
6(S  - 6S_)  = 0 

when  £ C 

(103) 

Physically  this  means  that  very  close  to  the  edges,  the  current  component 
parallel  to  the  edge  on  the  bottom  side  is  equal  to  half  that  on  the  top 
side  whereas  the  current  component  perpendicular  to  the  edge  on  the 
bottom  side  is  equal  to  zero. 
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Without  introducing  the  factor  5 ( S - 6S_)  which  is  known  as 
edge  modification,  we  cannot  interpret  the  unknown  J(r)  as  the  current 
density  on  the  top  side,  because  as  pointed  out  by  Jones  [16]  , the  inte- 
gral over  both  sides  of  the  plate  can  be  replaced  by  an  Integral  over 
one  side  provided  that  the  current  is  replaced  by  the  difference  between 
its  value  on  the  two  sides.  Therefore,  by  using  the  [l  — 6 ( S — 6S_)] 
factor  in  our  equation,  we  have  included  (although  not  rigorously)  the 
contribution  due  to  the  bottom  side  in  our  computation. 


4.7  Boundary  Conditions 


On  account  of  the  edge  condition,  we  require  that  the  normal 
component  of  the  current  must  vanish  on  the  rim  of  the  plate.  So  we 
impose  this  as  the  boundary  condition  in  our  solution.  Mathematically, 
we  have 


where 


A 

ne  *&(*)•  0 (104) 

ng  if  a unit  normal  to  the  edge  of  concern. 
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NUMERICAL  RESULTS  AND  DISCUSSIONS 

The  expressions  of  all  pertinent  equations  and  functions  as 
derived  in  Chapter  IV  are  programmed  on  a CDC  6500  system  at  TRW  Defense 
and  Space  Systems  Group.  A brief  description  about  the  computer  codes  is 
given  in  Appendix  B.  This  chapter  presents  the  numerical  results  obtained 
for  various  parameters. 

5.1  Flat-Square  Plate 

(a)  Current  Distributions  for  Various  Plate  Sizes 

(Less  than  One  Wavelength) 

Figures  8 and  9 show  the  plots  for  |jx|  along  the  x-  and  the 
y-axes  respectively  for  various  plate  sizes  at  normal  incidence  ($  - 0) . 

In  each  computer  run,  the  thickness  to  plate  length  ratio  was  kept  con- 
stant |V2oj  and  no  edge  conditions  were  imposed.  The  number  of  nodal 
points  used  is  16  with  four  nodes  in  each  x and  y direction.  It  is  seen 
that  in  general  the  shapes  of  |jxj  are  consistent  with  physical  intuition 
— that  is  for  small  plate  sizes,  |jx|  should  decrease  from  the  center  in 
the  x direction  and  Increase  in  the  y direction.  Since  the  same  number 
of  nodes  was  used  in  each  case,  it  is  logical  that  the  results  for  small 
plates  should  be  more  accurate  than  those  for  larger  plates.  The  plots 
for  the  phase  part  show  a fairly  constant  behavior.  Since  only  16  nodes 
were  used  in  each  direction,  the  current  does  not  show  a very  strong  edge 
behavior.  However,  Improvement  can  be  obtained  if  more  nodes  were  used 
as  shown  in  the  next  section. 

(b)  Convergence  Test 

The  numerical  solution  obtained  for  a well  posed  boundary-value 
problem  using  the  FEM  should  converge  to  the  exact  solution  as  the  ele- 
ment size  decreases.  This  is  illustrated  in  Figures  10  and  11  for  a 
square  plate  of  0.6X  in  size.  The  numerical  results  presented  here  were 
obtained  using  three  different  numbers  of  total  nodes  16,  25  and  36  (with 
equal  number  nodes  in  the  x and  the  y directions)  and  linear  shape 
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Figure  8.  Distributions  of  Jx  Along  x-Axls  for  Flat-Square  Plates  of  Various 
Dimensions  at  Normal  Incidence  ($  • 90°)  and  0-0. 
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Figure  9.  Distributions  of  Jx  Along  y-Axls  for  Fist-Square  Plates  of 
Various  Dimensions  at  Normal  Incidence  (+  • 90°)  and  • - 0. 
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Distributions  of  Jx  Along  x-Axls  for  s Flat-Square 
Plate  for  Various  Numbers  of  Nodal  Points,  Normal 
Points,  Normal  Incidence  ■ 90*)  and  9-0. 


I 
I 

function.  Here,  higher  Integration  accuracy  was  enforced  In  the  numeri- 
cal Integration.  Within  numerical  accuracy,  the  results  do  Indicate  con- 
vergence. However,  convergence  Is  not  very  uniform,  since  the  plots  for 
j Jx | show  more  disagreement  near  the  plate  center  than  elsewhere.  Careful 
examination  of  many  numerical  results  reveals  that  the  points  near  the 
plate  center  are  more  sensitive  to  Integration  accuracy,  plate  thickness, 
and  the  number  of  nodes  used.  One  reason  for  this  sensitivity  Is  due  to 
the  way  in  which  the  numerical  integration  is  performed.  In  a typical 
computer  run,  for  economical  reasons,  the  Integration  is  performed  over 
only  a specified  region  around  a given  node.  Therefore  the  nodes  close 
to  the  center  may  be  coupled  to  more  nodes  than  the  nodes  near  the 
edges. 

(c)  Plate  Thickness 

Figures  12  and  13  show  the  distributions  of  Jx  along  x-  and 
y-axes  respectively  for  a square  plate  of  0.4X  length  and  various  thick- 
ness (2t) . Sixteen  node  points  and  linear  shape  functions  were  used  for 
the  computation.  It  is  seen  that  the  amplitude  of  the  current  varies  in 
an  inverse  manner  to  the  plate  thickness.  As  the  plate  thickness  becomes 
larger  than  .02X,  some  small  ripples  and  a dip  occur  on  the  plots.  How- 
ever, the  phase  remains  relatively  constant.  As  will  be  discussed  in  a 
later  section,  the  plate  thickness  is  a relative  — not  absolute  — param- 
eter. It  is  found  that  the  numerical  solution  depends  actually  on  the 
ratio  of  the  plate  thickness  to  the  subdivision  size  in  the  integration. 

(d)  Linear  and  Sinusoidal  Shape  Functions 

Figures  14  and  15  show  the  comparison  of  |jx|  for  a square  plate 
of  0.6X  long  at  normal  Incidence  using  linear  and  sinusoidal  shape  func- 
tions. The  difference  between  these  two  results  is  small  for  the  size  of 
plate  considered  here.  But  for  larger  element  sizes,  the  difference  will 
be  more  noticeable,  and  it  is  expected  that  the  sinusoidal  shape  function 
produce  better  results.  However,  there  is  an  advantage  of  using  the  lin- 
ear shape  function  in  the  computation,  since  it  consumes  less  time. 


-.2  -.1  0 .1  *2 

y (X) 

• 

Figure  13.  Distributions  of  Jz  Along  y-Axis  for  Flat-Square  Plates  of 

Variouf  Plate  Thickness,  Normal  Incidence  (♦  “ 90°)  and  6 ■ 0. 
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(e)  Cross-Polarization  Current 


In  most  of  our  computation  the  polarization  of  the  Incident 
vector  was  assumed  In  the  x direction  (6*0);  hence  the  cross  polarization 
current  J is  much  smaller  than  the  parallel-polarization  component  J . 

y * 

Because  of  its  smallness,  it  is  difficult  to  obtain  reliable  results  with- 
out using  extremely  accurate  integration  scheme.  As  the  angle  of  polari- 
zation increases,  the  cross-polarization  current  will  likewise  increase. 
This  is  illustrated  in  Fig.  16  for  square  plate  of  0.4X  size. 

(f)  Boundary  Conditions 

It  is  known  that  in  solving  boundary-value  problems,  the  bound- 
ary conditions  must  be  included  in  order  to  obtain  a unique  solution. 
However,  in  variational  problems,  the  solution  automatically  satisfies 
the  natural  boundary  conditions.  In  our  study  here,  we  have  done  some 
investigation  about  the  effect  of  the  edge  condition.  The  edge  condition 
employed  here  is  given  by  n • J ■ 0,  where  n denotes  a unit  vector  nor- 
mal  to  the  edge.  Figure  17  shows  |JX|  for  a square  plate  of  0.7X  in 
length.  It  is  seen  that  the  difference  between  the  two  results  is  negli- 
gible in  the  y direction  and  significant  in  the  x direction.  This  dis- 
crepancy is  primarily  caused  by  the  few  number  of  nodes  used.  The  differ- 
ence will  probably  be  reduced  when  more  nodes  are  employed  instead  of 
only  16  being  used  here. 

(g)  Polarization  Angle 

Figures  18  and  19  show  the  effect  of  polarization  on  the  distri- 
bution of  Jx  along  x-  and  y-axes  for  a square  plate  of  0.6X  size  at  normal 
incidence.  As  expected  the  magnitude  of  Jx  is  proportional  to  cos  6 and 
the  phase  is  relatively  constant.  The  results  have  the  same  kind  of  behav- 
ior as  the  polarization  angle  on  the  cross-polarization  current  as  dis- 
cussed previously. 

(h)  Oblique  Incidence 

The  effect  due  to  angle  of  incidence  + on  the  current  distribu- 
tion is  illustrated  in  Fig.  20.  For  oblique  incidence,  we  have  symmetry 
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Figure  17.  Distributions  of  |J*|  Along  x-Axia  and  y-Axls  fcr  a Flat-Square  Plate 
With  and  Without  Edge  Condition,  Norael  Incidence  (♦  - 90°)  and  6-0. 
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Figure  18.  Distribution*  of  |Jx|  Along  x-Axia  for  s Flat-Square  Plate  for 
Various  Polarization  Angles  and  Normal  Incidence  ($  • 90°) . 
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only  with  respect  to  the  y-axis,  and  therefore  we  have  to  use  more  nodal 
points  in  the  y direction.  In  our  computation  here,  28  nodes  were  used 
with  4 in  the  x direction  and  7 in  the  y direction  for  each  oblique  inci- 
dence and  only  16  for  normal  Incidence.  It  is  seen  that  |JX|  along  x 
direction  increases  slightly  as  $ decreases.  The  results  for  normal 
incidence  and  oblique  incidence  differ  more  in  the  y direction,  especially 
in  region  close  to  the  edge.  The  difference  between  the  two  oblique  inci- 
dence cases  is  small.  From  the  classical  solution  for  the  diffraction  of 
a semi- infinite  plate  [l7],  we  know  that  the  current  in  the  neighborhood 
of  the  leading  edge  is  proportional  to  sin  £-jO-80°-$)J  . Our  numerical 
results  do  show  this  behavior. 

(i)  Current  Distribution  for  Resonant  Size 

Figures  21  and  22  show  the  current  distribution  of  Jx  in-  the  x 
and  y directions  respectively  for  a square  plate  of  one  wavelength  long 
and  various  thickness.  Even  though  only  16  nodal  points  were  used,  the 
numerical  results  are  surprisingly  good.  It  is  seen  that  the  results 
here  are  more  stable  with  respect  to  change  in  plate  thickness  than  the 
cases  for  small  plate  sizes.  The  value  of  |JX|  at  center  is  very  close 
to  the  physical  optics  prediction.  Again  the  phase  part  is  relatively 
constant  over  the  plate  except  near  the  edges.  The  distribution  is  found 
to  increase  a little  bit  faster  than  expected.  However,  if  more  nodal 
points  are  used,  the  distribution  of  |JX|  will  improve  as  shown  in  Fig.  23. 

5.2  Bent-Square  Plate 

Figure  24  illustrates  the  current  distribution  of  |JX|  along 
xc  direction  (y  • ,3X)  and  along  axis  for  a bent-square  plate  of 

- 130°  and  normal  incidence.  Again  it  shows  that  the  results  are  very 
close  for  the  two  different  shape  functions  used.  Comparing  these  results 
with  those  for  the  flat-square  plate,  we  find  that  the  bent  causes  changes 
on  the  current  distribution  mainly  in  the  region  close  to  the  bent.  The 
current  on  the  inclined  side  of  the  plate  is  larger  at  9 ■ 45° . This  is 
understandable  since  the  Inclined  side  is  at  normal  incidence  when 
9 - 45®.  The  distributions  of  |JX|  and  [ Jy(  along  xc  and  ^/y^  axes  are 
given  in  Figures  25  and  26  for  the  same  plate  at  normal  incidence  with 
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Figure  23.  Distributions  of  |JX|  Along  x-  and  y-Axes  for  a Flat-Square  Plate  with 
More  Numbers  of  Nodes  at  Normal  Incidence  ($  ■ 90°)  and  9 ■ 0. 
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Figure  24.  Distributions  of  jj^j  Along  x£  (y  - .3  X)  and  Along  y/yfe  Axis 
for  a Be.it-Square  Plate  ($.  - 130°)  for  Different  Angles  of 


Incidence 
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Figure  26.  Distributions  of  |Jy|  Along  and  y/yfc  Axes  for  a Bent-Square 
Plate  (<J»b  - 130°)  at  Normal  Incidence  (<j>  - 90°)  and  6-0. 
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polarization  angle  of  45°.  Unlike  the  flat-square  plate,  |jx|  here  is 
larger  than  |jy|.  This  indicates  that  the  presence  of  the  bent  reduces 
|jy|.  It  is  also  seen  that  the  derivative  of  |jy|  experiences  a discon- 
tinuity at  the  bent.  This  may  be  caused  by  the  existence  of  a line  charge 
there.  Figure  27  shows  the  distribution  of  |JX|  for  plate  size  of  one 
wavelength.  Two  angles  of  incidence  are  considered:  $ ■ 90°  and  45°. 

Again  as  for  oblique  incidence  28  nodes  were  used  in  each  computation. 

To  increase  matrix  stability,  a smaller  thickness,  A * 0.0375X,  was 
chosen.  It  is  seen  that  |JX|  in  the  x direction  is  similar  to  that  for  a 
flat-square  plate  while  |JX|  near  the  inclined  surface  in  the  ^/y  direc- 
tion becomes  smaller  and  shows  some  multiple  diffraction  phenomenon.  The 
value  of  |JX|  at  center  for  $ ■ 90°  is  very  close  to  the  physical  optics 
prediction. 

5.3  Computational  Time  and  other  Factors 

The  computer  time  required  in  each  run  depends  strongly  on  the 
number  of  nodes,  the  number  of  subdivisions  in  each  element  (three  sub- 
divisions were  used  in  most  cases  here) , the  integration  accuracy  as  well 
as  the  type  of  shape  function  used.  For  a 16  nodal  point  case,  the  CPU 
time  for  a typical  run  is  about  75  seconds.  Most  of  the  time  is  consumed 
in  matrix  fill.  Computational  time  increases  rapidly  with  the  number  of 
nodes.  In  order  to  reduce  the  computational  time,  it  is  proposed  that 
the  entire  double  surface  integral  be  transformed  into  a repeated  surface 
integral  in  the  variational  integral  equation.  The  transformation  is 
given  in  Appendix  A.  Another  idea  being  pursued  is  to  handle  the  singu- 
larity integration  analytically  rather  than  numerically  with  a finite 
thickness . 

Examination  of  the  numerical  results  shows  that  the  stability 
and  accuracy  of  the  numerical  solution  depend  strongly  on  many  parameters 
which  are  characterized  by  the  following  ratios:  (1)  thickness  (2t ) / 

subdivision  size  (6),  (2)  subdivision  size  (6) /subdivision  size  for  the 
singularity  Integration  (6S)  and  the  subdivision  size  (6)  itself.  Table  1 
shows  the  dependence  of  the  value  of  |jx|  at  the  plate  center  on  these 
parameters  for  a plate  size  of  0.6X  in  length. 


5-24 


*c  (*) 


H1  5- 


v / 

\ ,'! 

\ \ // 


SURFACE  DISTANCE  ALONG  y/yb  AXIS  (X) 
— y-AXIS  mim 


yb-AXis  - 


Figure  27.  Distributions  of  |JX|  Along  xc  and  y/y^  Axes  for  a Bent-Square 
Plate  ■ 130°)  at  Different  Angles  of  Incidence  and  6*0. 
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Table  1.  Dependence  of  |JX|  on  Plate  Thickness,  Subdivision 
Size  and  Element  Size  ($  - 90°,  0-0) 


A 

Element 
Size  (X) 

2t 

Thickness 

(X) 

6 

Subdivision 
Size  (X) 

6s 

Subdivision 
Size  for  Self 
Term  (X) 

lJ*l/ 

Vi 

Plate  Center 

.06 

.03 

.02 

.0075 

1.76 

.06 

.03 

.02 

.0075 

1.73 

.075 

.03 

.025 

.0125 

2.61 

.075 

.03 

.025 

.009375 

1.91 

.10 

.03 

.0333 

.0125 

2.34 

.10 

.03 

.0333 

.00833 

1.70 

.10 

.0225 

.0333 

.00833 

1.61 

.10 

.03 

.0333 

.0125 

2.34 

Experience  Indicates  that  to  obtain  good  results  we  should  choose  these 
parameters  carefully.  Unfortunately  there  are  no  exact  rules  for  us  to 
follow  in  choosing  these  parameters.  However,  It  is  recommended  that  in 
our  computation  we  should  observe  the  following  guidelines: 


* the  ratio  t/6  should  be  close  to  unity 

* the  subdivision  size  & should  not  be  larger  than  .031 

* the  ratio  */fi8  should  be  about  2.5. 
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CHAPTER  VI 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  problem  concerning  scattering  of  a plane  electromagnetic  wave 
from  a conducting  flat/bent- square  plate  has  been  solved  numerically  using 
the  finite  element  method.  The  singularity  of  the  Green's  function  is 
treated  by  assuming  a finite  thickness.  Two  computer  codes,  one  for  the 
flat-square  plate  and  the  other  for  the  bent-square  plate  were  written. 

Each  computer  code  can  handle  different  angles  of  incidence,  polarization 
angles,  and  plate  sizes.  Various  numerical  results  for  different  parameters 
were  presented  and  discussed.  It  is  found  that  for  plate  size  less  than 
or  equal  to  0.8  X,  the  numerical  results  are  convergent,  stable  and  good 
with  sixteen  nodes.  However,  when  the  plate  size  increases  beyond  this 
limit,  one  has  to  increase  the  number  of  nodal  points.  Another  unsatis- 
factory feature  is  the  thickness  of  the  plate.  To  circumvent  this,  a 
study  is  being  initiated  presently  to  treat  the  singularity  of  the  Green's 
function  analytically,  and  include  enough  consideration  for  the  current 
on  the  back  side  of  the  plate.  It  should  be  pointed  out  that  all  these 
problems  are  not  due  to  the  FEM  but  rather  it  is  inherent  in  the  E-field 
integral  equation  used.  This  suggests  that  all  these  difficulties  may 
be  avoided  if  other  more  suitable  equations  are  employed.  Investigation 
is  also  being  carried  out  in  this  area  of  alternate  formulation. 

Another  area  which  deserves  consideration  is  the  reduction  of 
the  double  surface  integral  into  a repeated  integral.  The  derivation  of 
this  has  been  given  in  Appendix  A.  This  approach  will  reduce  the  computer 
running  times  by  order  of  magnitude  as  has  been  demonstrated  by  TRW  for 
one-dimensional  problems.  Thus,  the  analytical  handling  of  the  singular- 
ity and  an  application  of  the  repeated  integral  technique  will  signifi- 
cantly improve  the  code  efficiency. 

The  FEM  method  is  completely  general  with  respect  to  geometry 
and  material  properties.  FEM  solutions  are  formulated  by  using  variational 
integrals  which  may  often  be  proportional  to  the  Lagrangian  or  Hamiltonian 
of  the  system  — thus  bearing  close  relation  to  the  physically  measurable 
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quantities.  The  technique  can  be  systematically  used  to  solve  such  com- 
plex and  difficult  problems  as  inhomogeneous  materials,  nonlinear  behav- 
ior, anisotropic  structures,  complicated  boundary  conditions  and  transient 
problems.  In  FEM,  the  body  is  discretized  without  destroying  either  the 
structural  details  or  the  continuous  nature.  No  separate  interpolation 
process  is  needed  to  extend  the  solution  to  every  point  within  the  contin- 
uum. Even  though  solution  is  obtained  at  a finite  number  of  discrete 
nodal  points,  the  formulation  inherently  provides  a solution  at  all  other 
locations  in  the  structure.  Boundary  conditions  are  very  naturally 
included  in  this  approach  and  only  the  geometrical  boundary  conditions 
need  be  specified.  The  method  gives  rise  to  a matrix  that  is  stable,  sym- 
metric and  often  banded. 


APPENDIX  A 


# 


TRANSFORMATION  OF  THE  DOUBLE  SURFACE  INTEGRAL  INTO  A 
REPEATED  SURFACE  INTEGRAL  IN  THE  VARIATION  INTEGRAL  EQUATION 


As  shown  in  Eq.  (34)  the  variational  integral  equation  involves 
a double  surface  integration.  From  a numerical  solution  point  of  view, 
this  is  a disadvantage  since  it  requires  large  computational  time.  One 
possible  way  of  reducing  the  computational  time  is  to  transform  the  double 
surface  Integral  into  a repeated  surface  integral  as  given  below. 


Let 


G(r|r')  • J(r')dS'  dS 


(A-l) 


J(r)  - Jx(x,y)x  + Jy(x,y)y  (A-2) 

From  (A-l),  (A-2)  and  using  the  coordinate  system  as  shown  in  Fig.  A-l,  we 
obtain 


AxxJx(x*y)Jx(x’’y,)  + XxyJx(x,y)Jy^x’,y^  + XxyJx^X' ,y' ^Jy(x,y) 

+ XyyJy(x',y')jy(x,y)J  dx'dy’dxdy  (A-3) 

X.  X and  X are  given  by  (37)  to  (39). 
xx  xy  yy 


Let  us  first  consider  the  first  integral  and  let 
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XX  X 
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n2b 

. 


2a.  2b 


y2y  \(x, y)  d,  dyjf  2 x„(x,,.«'  ,y’)  /,(«•  .y 


) dx’  dy' 


(A-4) 


• 2a  /*  2b 


Jx(x,y)dx  dy  n XX 


X (x-x'.y-y')  J-U’.y')  <**'  dy' 


Now 


let  us  change  the  variable  of  integration 


t - x-x 


u - y-y 

Substitution  of  (A-5)  into  (A-4)  yields 

2a  2b  x-2a  y-2b 


X (t,u)J  (x-t.y-u)dtdu 
XX  x 


Let 


♦JOC(t.«.*.y)  - x^ ( t . u) Jx (x-t , y-u) Jx (x » y ) 

By  changing  the  order  of  integration  (Figure  A-2)  we  obtain 


(A-5) 


(A- 6) 


(A- 7) 


2a+t  2b+u  -i 

.y;? ydu  | ' j j ♦xx(t»u»x»y)  dx  dyj 


XX 


-2a  -2b 
■2a/*  2b 


o o 


J^2y2b  ♦xx(t’u’x’y)  dX  dy  ] 


(A- 8) 
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If  we  let  t ■ -t  and  u * -u  In  the  first  integral  we  obtain 


2b  2a  2b 


2a- b 2b-u 


Ixx  mLir  \1  L (t,u,x,y)  */./.  (- 1 , -u , x , y )j  dxdy 


(A-9) 


Considering  the  integration  inside  the  bracket,  since  t and  u are  merely 
parameters  there,  we  can  let 

x ■ x + t , t**t  for  the  first  integral 

y ■ y + u , u ■ u for  the  second  integral 

We  get  from  (A-9) 

/-2a^2b  r /'2a-t/.2b-u  /*2a-t  /-2b-u  ~\ 

Ijn  mjdt J du  | J J ^(t.u.xft.y+u)  + J J ^(-t.-u.x.yjjdx  dy 

(A-10) 


Since  from  (A-7)  we  have 

^(-t.-u.x.y)  - X^(-t,-u)Jx(x+t,y+u)Jx(x,y) 
<t>  (t,u,xft,y+u)  - X (t,u)J  (x,y)J  (x+t,yfu) 


(A-ll) 


and  from  (37)  we  have 


Axx(t’u)  " Axx(_t*_u) 


- g(t.u)  + | y [(° + r)  - t]  -n  j 


1 e"3kr 

g(t,U)  - -fc-g “ 


and  R ■ J t2  + u2  + t 2 tfc  stands  for  the  half 

* thickness  now 


therefore  we  have 


$ (t,u,xft,y+u)  - ♦ (-t,-u,x,y)  and  (A-9) 

XX 
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becomes 
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In  a similar  manner  we  obtain 


r2a-t  r 2b-u 


Jx(xft,y+u)  Jx(*»y)  dx  dy 


•'o  o 


/*2a/*2b/~  2a/*2b  . . 

Iyy  “ / / / / ^yy  Jy(x.y)  Jy(x',y')  dx'  dy'  dx  dy 

./  o J o Jo  Jo 


a: 


rl&rlh  >.2b-u 

Tyy  “ 2 J q J Q Xyy(t>u)  dt  du  / / Jy(x*-t,y+u)  Jx(x,y)  dx  dy  (A-14) 


Next  let  us  consider  the  cross  product  terms 


rla.  r 2b/*  2a/-  2b  . . 

m / //  / Xxy  Jx(x,y)  Jy(*  »y  ) dx  dy'  dx  dy 

J o J o Jo  Jo 


(A-15) 


Following  the  same  procedure  as  before  we  obtain 


7*  2a  /2b  r r 2a-t  r 2b-u  /*  2a-t  r Zb-u  *| 

.*/.*■[/  /.  V*-"-’0 +/  J0 


2a-t  r 2b-u 


(A-16) 


However,  in  this  case  we  cannot  combine  the  two  integrals  inside  the  bracket 
together  because 


*xy(t,u,x+‘t,y+u^  " Xxy^t,u)  y*’y>  Jx(x+t,y+u) 
*(-t,-u,x,y)  - l(-t,-u)  Jx(x,y)  J (xft,y+u) 


(A-17) 
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sl"“  tu[f0  + s)2-f] 


So  4>  (**t,-u,x,y)  ” XTO(t,u)  J (x,y)  J (x+t,y+u) 

xy  *y  A y 

Substitution  of  (A-17)  and  (A-18)  into  (A-16)  we  get 


/•2a  * 2b  ( /*2a-t  >.2b-u 

mJ  J dtdu|y  y Xxy(t,u)  [ Jx(*ft»yfu)  Jy(x.y) 

+ Jx(x,y)  Jy(x+t,y+u)j  dx  dy | 


+ Jx(x,y)  Jy(x+t, 


Similarly  we  have 


/•  2a  /•  2b  ( /*2a-t  /•  2b—u  r 

lyx"/  / dt  du  < / / Xxy(t»u)  [ Jx(x+t,y+u)  Jy(x,y) 

•'O  •'O  ^ o o 


+ Jx(*»y)  Jy(x+t 


,y+u)J  dx  dy  j 


Combining  (A-13) , (A-14) , (A-19)  and  (A-20)  we  have 


n2b  /-2a-t  /•  2b-u 

A^(t,u)  dt  du  / / J,(x,y)  J^(x+t,y+u)  dx  dy 


/•za-t /-ZD-u 

1 1. 


-n 


2a  /-2b  , x , /*2a-t/-2b-u  , 

j Ayy(t,u)  dt  du  y j Jy(x,y)  Jy(»ft,y+u)  dx  dy 


+ 2X  TV(t’U)  dt  dU /2a  y 2b  U[J*(xft,3rfu)  Jy(x,y) 

+ Jy(x*t,yHi)  Jx(x,y)j  dx  dy  (A-21) 


which  is  the  repeated  surface  integral  form. 


since. 


From  a numerical  point  of  view  Eq.  (A-21)  is  superior  to  (A-3) 


1*  Uroit  of  integration  of  one  of  the  surface  Integrals 

is  reduced. 

2.  The  singularity  appears  only  in  one  of  the  integrals. 

3.  Since  in  the  repeated  surface  integral,  the  integrand 
involves  simple  functions,  such  as  polynomials  for  the 
linear  shape  function  and  sinusoidal  functions  for  the 
sinusoidal  shape  function,  the  integration  may  be 
integrated  out  analytically. 
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APPENDIX  B 


l 


A BRIEF  DESCRIPTION  OF  COMPUTER  PROGRAMS 

Since  a system  operators'  and  users'  manual  will  be  written 
later,  only  a brief  description  of  the  two  computer  programs  will  be  given 
here. 


1.0  Names  of  Programs 

Finite  (Benp-Square  Plate)  and  Element  (Flat-Square  Plate) 

2.0  Language 

FORTRAN  IV 

3.0  Input 3 

The  following  are  the  inputs  to  be  supplied  and  their  meaning: 
(a)  For  the  Flat-Plate  (PROGRAM  ''ELEMENT") 


DELT 

— 

thickness  of  plate  (X) 

AL2 

— 

half  length  of  plate  (X)  (y-axis) 

BL2 

— 

half  width  of  plate  (X)  (x-axis) 

KAX 

— 

number  of  the  divisions  for  the  width,  if  value  is 
less  than  2,  program  will  abort 

KAY 

— 

number  of  divisions  for  the  length,  if  value  is 
less  than  2,  program  will  abort 

JX 

— 

number  of  divisions  in  the  self-term  integration 
in  the  x direction 

JY 

number  of  divisions  in  the  self-term  integration 
in  the  y direction 

Typical  value  for  JX  and  JY  is  8 or  10. 

NKAX 

— 

number  of  subdivisions  for  the  singular  cell  inte- 
gration of  the  self-term  in  the  x direction 

NKAY 

number  of  subdivisions  for  the  singular  cell  inte- 
gration of  the  self-term  in  the  y direction 

Typical  value  for  NKAX  and  NKAY  is  10. 

JBC 

a control  parameter  for  boundary  condition.  If 
JBC  ■ 0,  no  edge  boundary  condition  is  imposed;  if 
JBC  ■ 1,  the  edge  boundary  condition  J • n "0  is 
used . 
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KSIN 

a control  parameter  for  choosing  shape  function. 
KSIN  = 0 for  linear  shape  function,  KSIN  = 1 
for  sinusoidal  shape  function. 

KMORE 

a control  parameter  for  integration.  KMORE  = 0 
for  normal  accuracy;  KMORE  ■ 1 for  additional 
accuracy. 

ZETA 

angle  of  polarization  of  the  E vector  with  the 
positive  x-axis,  0<  ZETA  < 90°. 

THETA 

— 

angle  of  incidence  made  with  the  negative  y-axis. 
0 < THETA  < 90° 

JCF 

a control  parameter.  JCF  «=  0 for  normal  inci- 
dence, JCF  = 1 for  oblique  incidence,  i.e., 
THETA  ^ 90°. 

NPRINT 

a control  parameter  for  printing  out  output.  If 
NPRINT  = 0,  normal  output  is  given;  NPRINT  = 1, 
more  output  such  as  the  matrix  is  dumped.  More 
output  is  needed  only  for  diagnostic  purpose. 

NDX 

■ 

number  of  subdivisions  for  each  nonsingular 
element  in  the  x direction. 

NDY 

— 

number  of  subdivisions  for  each  nonsingular 
element  in  the  y direction. 

Typical  value  for  NDX  and  NDY  is  3. 

XPOJ 

““ 

a control  parameter  for  integration 
Typical  value  is  2.0. 

XDEL 

a control  parameter  for  integration 
Typical  value  is  . 3X . 

For  the 

Bent-Plate  (PROGRAM  "FINITE") 

DELT 

— 

thickness  of  plate  (X) 

YL1 

— 

length  of  the  unbent  portion  of  plate  (X)  (y-axis) 

YL2 

— 

length  of  the  bent  portion  of  plate  (X) 

XL 

— 

half  width  of  plate  (X)  (x-axis) 

PHI 

— 

the  obtuse  angle  for  the  bent  90°  < <f>  < 180° , (deg) 

KAX 

— 

number  of  divisions  for  the  width,  if  value  is 
less  2,  program  will  abort 

KAY1 

— 

number  of  divisions  for  the  length  of  the  unbent 

portion,  if  value  is  less  than  2,  program  will 
abort 


JX,  JY,  NKAX,  NKAY,  JBC,  NPRINT , NDX,  NDY,  XPOJ, 
KSIN,  KMORE,  ZETA,  JCF,  THETA  and  XDEL  are  the 
same  as  described  in  the  flat-plate  case. 


4.0  Dimensions  Statement 
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In  using  the  programs,  all  the  arrays  should  be  dimensioned 
consistently  according  to  the  following  guidelines 

Let  L ■ number  of  elements  for  one  quadrant  of  plate 

M ■ total  number  unknown  quantities,  counting 
both  the  Jx  and  the  Jy  components. 

N = total  number  of  nodes 

J ■ total  number  of  subdivisions  for  nonsingular 
elements . 

We  have  for  the  flat-plate  case  — 

L = 2 (KAX)  (KAY) 

M - 2 (KAX+1)  (KAY+1) 

N - (KAX+1)  (KAY+1)  = M/2 

J =*  [(NDX)  (NDY)  + NDX]  /2  (assuming  NDX  «=  NDY) 
for  the  bent-plate  case  — 

L = 2 (KAX)  (KAY1  + KAY2) 

M = 2 (KAX+1)  (KAY1  + KAY 2 + 1) 

N - (KAX+1)  (KAY1  + KAY 2 + 1)  - M/2 

J = [(NDX)  (NDY)  + NDX]  /2 

and  the  dimensions  of  the  arrays  are  — 

A(M,M) , B(M) , W(M),  IZW(M),  CMAG (M) , PHS(M), 

ST (M,M)  , BS(M) , MLL(L),  MML(L) , MNL(L) , X(N), 

Y(N) , Z(N),  SUX(J,L),  SUY (J ,L) , SUZ(J,L), 

FSA(J,L) , FXY(NDX.NDY),  PMNE(N,L) 

(same  for  all  other  arrays  grouped  under  common 
block  C3) , CXX(L),  CYY(L),  CZZ(L),  TEX (NDX, NDY) , 
TEY(NDX.NDY),  and  TEZ (NDX, NDY) . 

5.0  Print  Out 

The  program  will  abort  if  the  value  of  KAX  or  KAY  in  PROGRAM 
ELEMENT  is  zero.  The  same  applies  to  Program  Finite  for  the  values  of 
KAX,  KAY1,  KAY2.  The  normal  print  outs  are:  the  supplier  input  values, 
the  node  numbers  and  nodal  current  densities  for  Jx  and  Jy  in  a|m2.  The 
current  densities  are  given  as  real  part,  imaginary  part,  magnitude  and 
phase . 
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6.0  Program  Flow 


The  main  program  performs  the  following  Items  in  order  — 

a.  Read  the  input  data. 

b.  Select  the  shape  functions. 

c.  Decide  whether  it  is  normal  or  oblique  incidence. 

d.  Calculate  the  x and  y coordinates  of  each  node 
and  other  essential  parameters  and  constants. 

e.  Call  in  subroutines  to  label  each  element  and  cal- 
culate the  area  coordinates.  All  these  values  are 
stored  in  common  blocks  for  later  uses. 

f.  Construct  the  linear  system  by  calling  in  sub- 
routines to  compute  the  singular  and  nonsingular 
matrix  elements.  To  save  computational  time,  the 
singular  matrix  element  is  computed  only  once. 

g.  Impose  the  edge  boundary  conditions  if  needed. 

h.  Invert  the  matrix. 

i.  Print  out  the  results.  A flow  diagram  is  shown  In 
Figure  B-l. 

7.0  Brief  Description  of  Subroutines 

XABEL  — This  subroutine  labels  the  triangular  elements  from 
t t 

one  to  N , N being  the  total  number  of  elements. 
Another  function  of  this  subroutine  is  to  identify 
the  correct  node  numbers  for  each  element.  Their 
values  are  stored  in  the  arrays  MLL,  MML  and  MNL  of 
the  common  block  TLAB.  It  is  called  only  once  in 
the  main  program. 

TRAN  — This  subroutine  computes  the  area  coordinate  param- 
eters, namely  P e,  Y e and  X 6 and  the  centroids 
mn  mn  mn 

for  each  element.  These  values  are  stored  in  two 
common  blocks,  C3  and  C18.  The  first  integer  in  an 
array  of  the  common  block  C3,  such  as  PMNE  represents 
the  node  number  while  the  second  integer  denotes  the 
element  number.  The  index  in  the  array  in  the  common 
block  C18,  such  as  CXX  represents  the  element  number. 
Another  function  of  TRAN  is  to  divide  each  element 


main  program 
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Figure  B-l.  Flow  Diagram  of  Programs 
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SUBROUTINE 

cor: 


SORT  OUT  THE 
* APPROPRIATE  AREA 
: COORDINATES  FOR 
-iTHE  ELEMENT  UNDER 
■ CONSIDERATION 


into  smaller  subdivisions,  either  rectangular  or  tri- 
angular in  shape.  The  number  of  the  subdivisions  is 
controlled  by  NDX  and  NDY  which  are  part  of  the  input 
data.  TRAN  then  computes  the  centroids  and  areas  of 
the  subdivisions  and  store  them  in  the  common  block 
TT.  TRAN  is  called  only  once  in  the  main  program. 

SEG  — This  subroutine  is  for  the  singular  element.  Because 
of  higher  accuracy  requirement  for  the  self-term  inte- 
gration (the  self-term  is  the  element  where  the  field 
and  the  source  points  coincide) , SEG  divides  the  sin- 
gular element  into  JX  and  JY  subdivisions  along  the 
x and  y directions.  Then  it  calls  in  subroutine  BAL 
to  do  the  computation.  To  save  computational  time, 
the  self-term  integration  is  only  performed  once  and 
the  data  are  stored  in  the  arrays  ST  and  BS  of  common 
block  SORT.  To  sort  out  the  correct  area  coordinate 
parameters  for  doing  the  computation,  SEG  calls  in 
CORD.  The  last  part  of  SEG  is  just  to  fill  up  the 
singular  matrix  element  for  each  node.  SEG  is  called 
only  once  in  the  main  program. 

BAL  — This  subroutine  computes  the  singular  matrix  element 
and  the  source  excitation.  The  source  excitation  is 
stored  in  the  vector  B.  To  calculate  the  dyadic 
Green's  function,  BAL  calls  subroutine  XKERN.  The 
results  from  XKERN  are  stored  in  the  common  block  AAA. 
BAL  then  uses  these  values  and  the  appropriate  area 
coordinates  as  represented  by  XCV,  XCW  and  XCQ,  etc. 
The  last  parameter,  KGO,  controls  computation.  If 
KGO  ■ 0,  computation  will  be  done  but  if  KGO  = 1, 
computation  will  be  excluded  and  previously  stored 
information  in  ST  and  BS  will  be  used  instead.  BAL 
is  called  only  in  SEG. 

CORD  — This  subroutine  has  a unique  simple  function  which  is 
to  sort  out  the  correct  area  coordinate  parameters 
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when  it  is  called  by  giving  the  node  and  the  ele- 
ment numbers  for  the  field  and  source  elements. 

The  resulting  information  is  stored  in  the  common 
block  C15.  This  common  block  has  eighteen  param- 
eters; nine  for  the  field  element  (the  ones  with 
P as  their  last  alohabet)  and  nine  for  the  source 
element.  CORD  is  called  both  in  SEG  and  TEG. 

STOR  — This  subroutine  is  called  in  the  main  program  when 
it  is  desired  to  compute  the  nonsingular  matrix 
elements.  The  primary  function  of  this  subroutine 
is  to  figure  out  the  number  of  elements  connected 
to  each  node.  From  Fig.  1,  it  is  understood  that 
the  number  of  elements  connected  to  a particular 
node  can  be  1,  2,  3 and  6 after  taking  symmetry 
into  account  (considering  only  one  quadrant  of  the 
plate) . STOR  then  calls  in  ADD  to  perform  the 
other  functions. 

ADD  — The  primary  function  of  this  subroutine  is  to  fig- 

ure out  how  the  field  element  and  the  source  ele- 
ment are  combined  as  far  as  a particular  node  is 
concerned.  Three  cases  can  arise:  (a)  both  the 

field  and  the  source  elements  are  connected  to 
that  node;  (b)  only  the  field  element  is  connected 
to  that  node;  in  this  case,  the  control  parameters 
are  set  by  letting  MSON  ■ 1 and  MDUB  =0;  (c)  only 
the  source  element  is  connected  to  that  node;  here 
we  have  MSON  ■ 0 and  MDUB  * 1.  For  each  combina- 
tion of  the  source  and  the  field  elements,  ADD 
then  calls  in  TEG  to  do  the  Integration.  ADD  is 
called  only  in  STOR. 

TEG  — This  is  a very  simple  subroutine.  Its  function  is 

to  divide  each  source  element  into  NDIM  subdivi- 
sions and  call  in  CAL  to  do  the  integration. 
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— For  the  nonsingular  matrix  elements,  this  subrou- 
tine does  the  most  important  job.  It  is  called 
TEG.  The  subroutine  makes  use  of  the  symmetry  of 
the  problem  by  dividing  the  area  of  the  plate  into 
four  quadrants.  To  save  computational  time,  accu- 
racy of  integration  is  controlled  by  several  input 
parameters  such  as  XPOJ,  XDEL,  NDX  and  NDY.  To  do 
the  computation,  CAL  calls  in  XKERN  to  do  the 
dyadic  Green's  function  computation  and  PART  to  do 
the  area  coordinates  computation  and  filling  in 
the  matrix  elements.  The  resulting  information 
from  PART  is  stored  in  the  common  block  SAD  which 
is  then  used  to  fill  in  the  main  matrix  A. 

— The  primary  function  of  this  subroutine  is  to 
handle  the  dyadic  Green's  function.  It  is  called 
in  BAL  and  CAL.  The  computed  values  are  stored  in 
the  common  block  AAA. 

— This  subroutine  is  called  only  in  CAL  and  its  func- 
tion is  to  calculate  the  area  coordinate  parameters 
and  their  combinations  with  the  dyadic  Green's 
function  as  computed  in  XKERN.  The  resulting  com- 
puted values  are  stored  in  the  common  block  SAD. 

In  computing  the  area  coordinate  parameters,  atten- 
tion is  given  to  the  way  the  source  element  and  the 
field  element  are  combined.  This  is  controlled  by 
the  control  parameters  MSON  and  MDUB  of  the  common 
block  TCH. 

— This  subroutine  is  sometimes  needed  to  do  more 
accurate  integration  by  subdividing  each  subdivi- 
sion again  into  smaller  subdivisions.  Because  of 
the  computational  time  involved,  this  subroutine  is 
not  used  most  of  the  time  for  smaller  plate  sizes. 
It  is  called  only  in  CAL. 
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This  subroutine  determines  whether  a given  integer 
is  even  or  odd.  It  is  called  in  XMORE. 


CSIMEQ  — This  is  a matrix  inversion  subroutine.  It  simply 
performs  the  inversion  of  a given  complex  matrix 
by  Gaussian  elimination  alogrithm.  In  calling  the 
subroutine,  the  first  integer  number  must  be  equal 
to  the  dimension  of  the  matrix  A;  the  second  inte- 
ger represents  the  number  of  unknowns  of  the  linear 
system;  the  third  parameter  is  the  matrix  to  be 
inverted;  the  fourth  parameter  is  a vector  to  store 
the  solution  after  inversion;  the  fifth  parameter 
is  the  excitation  vector  and  the  last  parameter  is 
a dummy  vector  for  storage  purpose. 

CABS1  — Its  function  is  to  determine  the  absolute  value  of 
a complex  number.  It  is  called  by  CSIMEQ, 

XSIN  — This  subroutine  is  to  compute  the  sinusoidal  shape 
function.  In  computing  the  sinusoidal  coordinates, 
the  subroutine  makes  use  of  the  area  coordinates 
stored  in  common  block  C15.  It  is  called  in  ADD, 
BAL  and  PART.  In  calling,  the  cartesian  coordin- 
ates of  the  observation  and  field  points  are  sup- 
plied and  the  resulting  values  are  returned  through 
arguments. 


8.0  System  Requirement 

A FORTRAN  compiler  w/standard  input/output  and  mathematical 
macros  is  all  that  is  needed  to  run  this  job.  The  storage  requirement  for 
running  a particular  job  depends  on  the  dimension  size  set  by  the  program- 
mer. For  a particular  problem  run  here  at  TRW,  it  took  90  cpu  seconds  and 
50  k of  core  to  run  and  execute  on  a CDC  174  computer. 
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